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A CONSTRUCTIVE ARBITRARY-DEGREE KRONECKER 
PRODUCT DECOMPOSITION OF TENSORS 

KIM BATSELIER AND NGAI WONG* 


Abstract. We propose the tensor Kronecker product singular value decomposition (TKPSVD) 
that decomposes a real k -way tensor A into a linear combination of tensor Kronecker products 
with an arbitrary number of d factors A = SjLi a j 0 • • • 0 A^ . We generalize the matrix 

Kronecker product to tensors such that each factor Aj J in the TKPSVD is a k -way tensor. The 
algorithm relies on reshaping and permuting the original tensor into a d-way tensor, after which 
a polyadic decomposition with orthogonal rank-1 terms is computed. We prove that for many 
different structured tensors, the Kronecker product factors A^\ ..., A^ are guaranteed to inherit 
this structure. In addition, we introduce the new notion of general symmetric tensors, which includes 
many different structures such as symmetric, persymmetric, centrosymmetric, Toeplitz and Hankel 
tensors. 


Key words. Kronecker product, structured tensors, tensor decomposition, TTrlSVD, general¬ 
ized symmetric tensors, Toeplitz tensor, Hankel tensor 
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1. Introduction. Consider the singular value decomposition (SVD) of the fol¬ 
lowing 16 x 9 matrix A 
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into a 16 x 16 orthogonal matrix U, a 16 x 9 diagonal matrix S and a 9 x 9 orthogonal 
matrix V. The distinct entries of A were drawn from a standard normal distribution. 
Note that A has only 5 distinct columns and 7 distinct rows, limiting the rank to 
maximally 5, and has no apparent structure such as a symmetry or displacement 
structure (Hankel or Toeplitz matrix). Reshaping the first left and right singular 
vectors into a 4 x 4 and 3x3 matrix respectively results in the following Hankel 
matrices 
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In fact, it turns out that all singular vectors of A can be reshaped into square Hankel 
matrices. In the same vein, other matrices can be easily constructed such that their 
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reshaped singular vectors are square structured matrices (symmetric, Toeplitz, cen- 
trosymmetric, persymmetric,...). And as we will demonstrate in this article, the same 
can even be done for k- way tensors. Now what is going on here? Why are all singular 
vectors so highly structured? This is explained in this article, where we introduce the 
tensor-based Kronecker product (KP) singular value decomposition (TKPSVD). The 
TKPSVD decomposes an arbitrary real k- way tensor A G R raiX " 2 x---xn k as 

R 

(1.2) A = 

3 =i 

where ® denotes the tensor Kronecker product, defined in Section [3j and the tensors 

Af G R"i i)x - x 4 i) satisfy 

d 

(1.3) = 1 with 4’ 1 = n r (r G {1,..., k}). 

i-1 


It turns out that the TKPSVD is key in explaining the structured singular vectors of 
the matrix A. In Section [lj we show that computing the SVD of (1.1) is an interme¬ 
diate step in the computation of the TKPSVD of a 12 x 12 Hankel matrix into the 
Kronecker product of a 4 x 4 with a 3 x 3 matrix. 


The Kronecker rank (5| , which is generally hard to compute, is defined as the minimal 
R required in (|1.2[) in order for the equality to hold. If for any r G {1, ...,&} we have 


that n = = rif 1 = 


= ni d) 


then n r — n d . For this reason, we call the number 


of factors d in (1.2) the degree of the decomposition. The user of the TKPSVD al- 


(i) 

gorithm is completely free to choose the degree d and the dimensions n r of each of 


the KP factors, as long as they satisfy (1.3). In addition to the development of the 
TKPSVD algorithm, another major contribution of this article, shown for the first 
time in the literature, is the proof that by our proposed algorithm we will have the 
following favorable structure-preserving properties when all A^ and A are cubical: 


if A is 
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The fact that each of the factors A';- 1 inherits the structure of A is not trivial. In 
providing this proof a very natural generalization of symmetric tensors is introduced, 
which we name general symmetric tensors. In addition, Toeplitz and Hankel tensors 
are also generalized into what we call shifted-index structures , which are special cases 
of general symmetries. In fact, when A is general symmetric, then all its cubical factors 
A 1 ;- 1 will also be general symmetric. Another interesting feature of the TKPSVD 


algorithm is that if the summation in (1.2) is limited to the first r terms, then the 


relative approximation error in the Frobenius norm is given by 


(1.4) 
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Equation (1.41 has the computational advantage that the relative approximation er¬ 
ror can be easily obtained from the aj’s without having to explicitly construct the 
approximant. 


1.1. Context. The TKPSVD is directly inspired by the work of Van Loan and 
Pitsianis [TQJ. In their paper, they solve the problem of finding matrices B,C such 
that \ \A — B ® C\\p is minimized. The globally minimizing matrices B , C turn out to 
be the singular vectors corresponding with the largest singular value of a particular 
permutation of A. In [5], the full SVD of the permuted A is considered and the 
corresponding decomposition of A into a linear combination of Kronecker products is 
called the Kronecker product singular value decomposition (KPSVD). Applications of 
the d = 2, k = 2 KPSVD approximation in image restoration are described in mm , 
whereas extensions to the d = 3, k = 2 case using the higher order singular value 
decomposition (HOSVD), also for imaging, are described in pT] [IT]. 

The decomposition in this paper is a direct generalization of the KPSVD to an 
arbitrary number of KP factors d and to arbitrary k- way tensors. In the TKPSVD 
case, the SVD is replaced by either a canonical polyadic decomposition (CPD) [21 
[6] with orthogonal factor matrices, the HOSVD [4] or the tensor-train rank-1 SVD 
(TTrlSVD) [Tj. The TKPSVD reduces to the KPSVD for the case d = 2 and k = 2. 
Contrary to previous work in the literature, we are not interested in minimizing 

u-T,U< 




Instead, we are interested in a full decomposition 


such that any structure of A is also present in the KP factors. Explicit decomposition 
algorithms for when d > 4 and k > 3 are not found in the literature. To this end, 
our proposed TKPSVD algorithm readily works for any degree d and any tensor 
order fc, does not require any a priori knowledge of the number of KP terms, and 
preserves general symmetry in the KP factors A^\ Furthermore, a Matlab/Octave 
implementation that uses the TTrlSVD and works for any arbitrary degree d and 
order k can be freely downloaded from https://github.com/kbatseli/TKPSVD, In 
brief, the contributions of this article are: 

• an explicit formulation of the generalization of the KPSVD algorithm for 
d > 2 and k > 2 is presented for the first time in the literature, 

• a new notion of general symmetric tensors, which describes many tensor struc¬ 
tures, is introduced, 


• we prove that for a general symmetric tensor A all KP factors A':^ in (1.2) are 
guaranteed to have the same general symmetry under a particular condition. 

The outline of this article is as follows. In Section [2] we introduce some basic 
tensor concepts and notations. In Section [3] we generalize the matrix Kronecker 
product into the tensor Kronecker product and present some of its properties. In 
Section [4] we first derive the theorem that underlies the TKPSVD and then present 
the TKPSVD algorithm for both general and diagonal tensors. In Section [5] we 
introduce the framework of general symmetry that describes many different structured 
tensors such as symmetry, centrosymmetry, Hankel, Toeplitz, etc.. Preservation of 
general symmetry in the KP factors when the original tensor is general symmetric is 
proven in Section [bj Numerical experiments that demonstrate different aspects of the 
TKPSVD are presented in Section [7] after which conclusions follow. 


2. Tensor basics and notation. We only consider real matrices and tensors 
in this paper. Scalars are denoted by greek letters (a,/3,...), vectors by lowercase 
letters (a, b ,...), matrices by uppercase letters ( A,B ,...) and higher-order tensors by 
uppercase caligraphic letters (A,B,...). The notation (-) T denotes the transpose of 
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either a vector or a matrix. A /cth-order or k -way tensor is a multi-way array A £ 
jgmixn 2 x---xn fc t ensor j s cubical if all dimensions are equal, e.g., ni = **2 = • • • = ***,. 

Entries of tensors are always denoted with square brackets around the indices. 
This enables an easy way of representing the grouping of indices. Suppose for example 
that A is a 4-way tensor with entries -4.[?;,][i 2 ][?: 3 ][i 4 ]- To improve readability we do not 
write the square brackets when all indices are considered separate, therefore A t1 ,; 2! ; 3! ; 4 = 
A[ii][i 2 ][i 3 ][i 4 ]. A 3-way tensor can now be formed by grouping, for example, the first 
two indices together. The entries of this 3-way tensor are then denoted by -A[ 2l 22 ][i; 3 ] [?, 4 ], 
where the grouped index [ 11 * 2 ] is easily converted into the linear index *1 + **i (*2 — 1 )- 
Grouping the indices into the [i\] and [ 12 * 3 * 4 ] results in a r*i x ** 2 * 13*14 matrix with 
entries -4[i 1 ][i 2 i 3 i 4 ]- The column index [* 2 * 3 * 4 ] is equivalent to the linear index *2 + 
** 2(*3 ^ 1) +** 2 ** 3(*4 — 1)- A very special case of grouping indices is obtained when all 
indices are grouped together. The resulting vector is then called the vectorization of 
A, denoted vec(A), with entries -A[i 1 i 2 i 3 i 4 ]- 

The r-mode product of a tensor A £ K n i x " 2 x -x«i= with a matrix U £ R pxrar is 
defined by 


(A x r U)i 1 ...i k _ 1 ji k+1 ...i d 

so that Ax r U £ R"ix-xnr-ixrxr. r+ ix-x nji The multiplication of a k-way tensor A 
along all its modes with matrices P \,..., P*. 

B = A xi Pi x 2 P 2 x 3 • • • x fc P fc , 


— ^ ^ j i, A ?. 1 ■ ■ ■ /' ,■ 


can be rewritten as the following linear system 

( 2 . 1 ) vec (B) = (Pd ® ■ • • <g> P 2 ® Pi) vec(A), 

where ® is the conventional matrix Kronecker product m, which is defined and 
generalized to the tensor case in Section [3] The inner product between two tensors 
A,B £ R niX -x"J is defined as 

(A,B) = ^2 Bhi 2 -id = vec(A) T vec(B). 

*1 )■“ )*eZ 


Two tensors A, B are orthogonal with respect to one another when (A, B) = 0. The 
norm of a tensor is taken to be the Frobenius norm ||A||f = (A, A) 1 ^ 2 . A k -way 
rank -1 tensor A £ R«i x '" x «fc i s per definition the outer product [ 8 ], denoted by o, of 
k vectors a W £ (* £ { 1 ,..., k}) such that 


( 2 . 2 ) 


A = o a^ 


o a 


(fe) 


with Ai 


(i) ( 2 ) 

= a h a l 


,(fc) 


Any real k -way tensor A can always be written as a linear combination of rank -1 
terms 


(2.3) 


A = J2a j afKafK. 
l=i 


o a 


(k) 


where < 7 j £ K and all the vectors satisfy ||a ^||2 = 
called a polyadic decomposition (PD) of the tensor A. 


1. Such a decomposition is 
When the equality in (2.31 
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holds with a minimal number of terms R , then the PD is called canonical (CPD) [2][6]. 
The number R in the CPD is called the tensor rank. Unlike the SVD, each of the 
rank-1 terms in the (C)PD is not necessarily orthogonal. 

The Tucker decomposition m |18] writes a k -way tensor A as the following 
multilinear transformation of a core tensor S £ M ri xr 2 x---xr k j-,y f ac ^ or matrices 


(2.4) 


A = S X! £/ (1) x 2 U (2) x 3 ••• x fe U (k \ 


which can also be written as (2.3) where each aj is now an entry of the core tensor 


S. Each rank-1 term of the Tucker decomposition is then given by 
S hi2 ... ik UW (:,<!) o tf< 2 >(:,* 2 ) o • • • o C/W(:,4), 
where we use Matlab colon notation to denote columns of the factor matrices. 


The minimal size of the core tensor S such that the equality in (2.4) holds is called 


the multilinear rank. The higher-order SVD (HOSVD) ^ is a Tucker decomposition 
where the core tensor S has the same dimensions as the original tensor A and with the 
additional property that the factor matrices U <J ' 1 and the slices of S in the same mode 
are orthogonal. This implies that each rank-1 term is orthogonal to all other rank-1 
terms in the HOSVD, which has the immediate advantage that the approximation 


error can be determined as in (1.4). 


The PARATREE/TTrlSVD decomposition [TJ 13] is another decomposition of a 
A:-way tensor into orthogonal rank-1 terms as described by (2.3). The total number 
of terms in the TTrlSVD is upperbounded by R = Tlr=0 min (n r +i, rii=r+2 ni) and 
therefore depends on the ordering of the indices. This decomposition is computed 
from repeated reshapings and SVD computations and is unique for a fixed order of 
indices. Note that although each of the rank-1 terms is orthogonal with respect to all 
others, unlike the HOSVD, the factor matrices U^ are not orthogonal. In addition, 
the scalar aj coefficients obtained in the PARATREE/TTR1SVD decomposition are 
guaranteed to be nonnegative. This has the advantage that one can plot these aj 
coefficients in descending order and inspect the relative weight of each of the rank-1 
terms in a very straightforward manner, just like one can do with the singular values 
of a matrix. 

3. Tensor Kronecker product. 

3.1. Definition. The definition of the Kronecker product for two matrices is 
well-known. Ii B £ R m i xm 2 anc j (j g jjruxn 2) then their Kronecker product B ® C is 
an mi x m 2 block matrix whose (i 3 ,* 4 )th block is the ni x n 2 matrix B l[t . l4 C 


/ 

( Bn ■ 

B\ ni \ 

1 1 

( BnC ■■ 

■ B lni C\ 

(3.1) | 


B ) 


{■B mil C ■■ 

Bmini CJ 


Generalizing this definition to the Kronecker product of k -way tensors is quite straight¬ 


forward, although not in the form as given by (3.1). Before giving the definition of the 


tensor Kronecker product, we first investigate how the entries of the matrix Kronecker 
product are described by the indices of the original matrices. The following lemma is 
easily verified. 
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Lemma 3.1. If the entries of the matrices B £ j£ TO ixm 2) (j g |"i x " 2 are denoted, 
B isii and C^i 2 respectively , then the entries of their Kronecker product A = B ® C 
are described by A[ ili3 ][ i2ii ] = Bi 3 i i Ci x i 2 for all possible values of i\, ii, 73,74. 

Remember that the grouped indices [*1*3], [ 12 * 4 ] are easily converted into the linear 
row index i\ + 711(73 — 1) and the linear column index i 2 + 712(74 — 1) respectively. The 
definition of the tensor Kronecker product follows from the generalization of Lemma 
|3.1| to multiple indices. 

Definition 3.2. Let B £ Rnixiux-xn^ e ^m 1 xm 2 x--xm k tw0 ^_ way ^ en _ 

sors with entries denoted by Bi k+1 ...i 2k ,Ci 1 ...i k respectively. The tensor Kronecker prod¬ 
uct A = B ®C £ R n i m iX“ 2 ™ 2 x-xni ! m t defined from 

■'4[iiifc + i][i 2 ifc+2]-” [iki2k] 

which needs to hold for all possible values of i \,..., i 2 k- 

Although the tensor Kronecker product is a very straightforward generalization 
of the matrix Kronecker product, we failed to find any reference to it in the litera¬ 
ture. One possible implementation of the tensor Kronecker product would be to use 
Definition |3.2| over all possible values of the indices i \,... ,frk but this would not be 
very efficient. Instead, one can use an existing implementation of the vector/matrix 
Kronecker product (‘kron.rn’ in Matlab) on the vectorized tensors vec(P), vec(C). In¬ 
deed, the entries of c = vec (B) ® vec (C) are indexed by the single grouped index 
[ii ■ ■ ■ ikik+i • • • * 2 &]. One can then reshape and permute the entries of c such that the 
desired [iiik+i] [* 2 *fc+ 2 ] ■ • • [ikfrk] index structure is obtained. This is how the tensor 
Kronecker product is implemented in our Matlab/Octave TKPSVD package. 

3.2. Properties of the tensor Kronecker product. We briefly list some 
properties of the tensor Kronecker product without going into details. The following 
properties are easily verified 

A®(B + C)=A®B + A®B, 

(A + B)®C = A®C + B®C, 

( 1 aA ) ® B = A ® (aB) = a(A ®> B), 

(A®B)®C = A®{B®C), 

where A, B , C are k -way tensors and a is a scalar. Just as in the matrix case, the 
tensor Kronecker product is not commutative but permutation equivalent. That is, 
there exists permutation matrices Pi,..., Pk such that 


(3.2) 


A®B = ( B®A ) xi Pi x 2 * • • x fc P fc . 


This is easily seen from the definition. Indeed, suppose we have that C = A® B and 
C = B ® A, then 


C[iii J . + i][i 2 ifc + 2]---[iki2fc] •Ai k+1 ...i 2k Bi 1 .„i k , 

C[i fc+ iii][i fc+2 i2]-"[i2feifc] = •Ai 1 ...i k Bi k+1 ...i 2k . 


Now let Pi,...,Pfc be the permutation matrices that swap [ik+iii] into [iifr+i], 
[*fc+ 2 * 2 ] into [* 27 / 0 + 2 ]) ■••) [frufr] into [ik^k] respectively, then (3.21 follows. Further¬ 
more, if A , B are cubical and of the same dimension then C and C are permutation 
similar, which means that P\ = Pz = ■ ■ ■ = P/o- 
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The mixed-product property of the Kronecker product states that if A, B , C, D are 
matrices such that one can form the matrix products AC, BD, then (A®B)(C(l)D) = 
(AC)®(BD). This is called the mixed-product property, because it mixes the ordinary 
matrix product and the Kronecker product. We can also write the mixed-product 
property using the 1-mode product as (C®D) Xi(A®B) = (Cxid)®(DxiB). Its 
generalization involves the mixing of the r-rnode product with the tensor Kronecker 
product. Let A, B be matrices and C,V be k- way tensors of appropriate dimensions 
then for any r £ {1,..., k} 

(3.3) [C. <g> V) x r {A®B) = (C x r A) (g> (V x r B). 


What the mixed-product property tells us is that we can obtain the r-mode product 
of the tensor C ®T> with the matrix A® B from the Kronecker product of (C x r A) 


with (V x r B ). Choosing r = 1 and replacing C,V with matrices in (3.3) results in 
the matrix mixed-product property. 


4. TKPSVD Theorem and algorithm. Using Definition |3.2| we can easily 
extend the tensor Kronecker product to multiple factors. Suppose we have three 3- 
way tensors A^\ A^, A^, then their Kronecker product A = <g> A ^ <g> A^ is 

completely characterized by the following relationship 


(4.1) 


Ar = 4« A (2) A (3) 


2122^3 ' 


Now suppose we permute all entries such that -A[ij* 2 i 3 ][* 4 * 5 * a ][* 7 * 8 * 9 ] = 


and that this is a rank-1 tensor. According to (2.2), A can then be written as the 
following outer product of vectors A = o a ^ o a ^ with 


(4-2) •' 4 [* i * 2 * 3 ][* 4 * 5 * 6 ][* 7 * 8 * 9 ] a ’- .---1 a ’~ - - 1 a <- 


(1) a (2) a (3) 

[* 1 * 2 * 3 ] [* 4 * 5 * 6 ] [* 7 * 8 * 9 ]' 


Comparison of (4.1) with (4.2) allows us to conclude that a*- 1 - 1 = vec(A^ 1 - ) ),= 
vec(A ( ' 2 )) and a®" = vec(A^ 3 ). We formalize this observation in the following Theo¬ 
rem. 

Theorem 4.1. For a given k-way tensor A £ f n i xn 2 >< -x ,i )= j if 


R 


(4.3) A — aj A 


(d) 




( 2 ) 


)A 


^ and A = aj o o • • • o a\ d \ 


l=i 


1=1 


where A is the permutation of A such that the indices of the vectors are identical 

to those of the k-way Aj 
{1,... ,7?}. 


tensors, then = vec(A'f > ) for all i £ {1,..., d},j £ 


(*)\ 


Observe that the order of the Kronecker products in (4.3) is reversed with respect 


to the order of the outer products. Theorem |4.1| is c rucial for the TKPSVD algorithm, 
since it tells us that the desired decomposition (1.21) can be computed from a PD of A. 


We now derive the TKPSVD algorithm by means of a simple example. Suppose we 
have a 3-way tensor A for which we want to find a degree-3 decomposition. This im¬ 


plies, as shown in (4.1), that each entry of A is labeled as -4.[i 1 * 4 ,i 7 ][* 2 i s * a ][* 3 *6* 9 ]• Figure 
|1. 1 ['a) illustrates how the grouped indices of the tensor A relate to those of the KP 


factors A, . Theorem 


4.1 


tells us that the desired TKSPVD can be obtained from a 


PD of the permuted tensor A. The first step in the TKPSVD algorithm is then to per¬ 
mute the indices of A such that their order corresponds with U, * 2 , * 3 , H, * 5 , H, * 7 , *8, * 9 - 
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4^71 4L 
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( 3 ) 




( 2 ) 


X, 


( 1 ) 


(a) = 3, d = 3 


If 




i. 

I 










14 i r * 2 

A" 







I 








=z-4M®JM®W 


(b ) k = 2, d = 2 (c) k = 2, d = 3 

Fig. 4.1: How the grouped indices of A relate to the indices of the KP factors Aj. 


In order to do this, we first reshape the 3-way tensor A into the 9-way tensor A with 
entries •/ 4 i 1 * 4 i 7 i 2 tBi 8 » 3 » 6 *g‘ The indices of A are then permuted into the desired order 
•' 4 iii 2 * 3 t 4 t 5 * 6 i 7 * 8 * 9 - The next step of the TKPSVD algorithm is to compute the KP 
factors A { j ] . each of which is computed as a vector in a PD. We therefore group the 
indices such that we obtain the 3-way tensor A with entries ^4[z 1 i 2 i 3 ] [z 4 i 5 i e ] [i 7 i s i 9 ] ■ The 
steps prior to the computation of the PD are hence summarized as 


A 


[212427] [22252s] [232629] 


reshape . 

^■‘^1242722^5^8^3^6^9 


permute 


>A 


2l2223242526«7*8 4 9 


reshape 


>A 


[212223] [2425261(27^8^9] ’ 

Each of the KP factors A!:- 1 is obtained from reshaping the a - ■* vectors of the PD (Fig¬ 
ure 4.2) into a 3-way tensor of the correct dimensions. In order to make sure this 


procedure of reshaping and applying the permutation is clear, we also demonstrate it 
for a simple matrix example. Suppose we have a 12 x 12 Hankel matrix A, which we 
want to decompose into a sum of KPs of a 4 x 4 matrix with a 3 x 3 matrix. If the 
entries of the KP factors A^,A^ are labeled by respectively, then the row 

index of A is [ 11 * 3 ] and the column index is [* 2 * 4 ], shown in Figure |4T| (b). The steps 
prior to the computation of the PD are now 


A 


2123 ] [ 2224 ] 


reshape 


^■^212 


permute 


y >Ai 


reshape 


> A\ 


t" 12'2 ] [23 2 4 ] ' 


The dimensions of the tensors in each of these steps are 12x12,4x3x4x3,4x4x3x3 
and 16x9 respectively. The 16x9 matrix A shown on the front page is in fact obtained 
from this procedure. The Hankel structure is apparently lost due to the consecutive 
reshapings and permutation. The final step is to compute a PD with orthogonal rank- 
1 terms, which for the matrix case is the SVD. The previous two examples might give 
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[* 7 * 8 * 9 ] [* 7 * 8 * 9 ] 




Fig. 4.2: Decomposition of the 3-way tensor A into a linear combination of rank-1 
terms. 


the impression that d and k need to be equal in the TKPSVD. This is not the case. 
In Figure [4~i~[ c) we show a degree-3 TKPSVD of a matrix. Here, the steps prior to 
the computation of the PD are 


A 


[* 1 * 3 * 5 ] [* 2 * 4 * 6 ] 


reshape 




permute 


21*3*5*22426 


> A, 


reshape 


21*2232425*6 


> A 


[* 1 * 2 ][* 3 * 4 ][* 5 ® 6 ]‘ 


The pseudo-code for our general TKPSVD algorithm is presented in Algorithm |4.1| As 
we will show in Section [6] the structure-preserving property of the TKPSVD critically 
depends on the fact that the rank-1 terms of the computed PD are orthogonal with 
respect to one another. Another consequence of this orthogonality is that 


11-^1 \f — \J H + erfj, 

where R is the total number of terms in the decomposition. The relative approxi¬ 
mation error obtained from truncating the number of terms to r < R is then easily 
computed as 


\\A ~ Y!j= 1 a 3 ^ ® ® A^Wf _ \/ a r+l “I I" a R 

II^IIf V a i^ - a R 

This is especially convenient when using the TTrlSVD [Tj to compute the PD, since 
then all ays are positive and can be sorted in descending order just like singular values 
in the matrix case. 


Algorithm 4.1. TKPSVD Algorithm 

Input: tensor A, dimensions n ^,..., nf '*, n ^,..., n'A ,..., n 

Output: cri,..., gr, tensors A^ 

A 4— reshape A into a (kd)-way tensor according to n{ \ ..., nf\ n A \ ..., 

A 4— permute A to indices 21 * 2 * 3*4 • • • ikd-Akd 

A 4— reshape A into a d-way tensor by grouping every k indices together 

..., \ cri,..., ctr 4 — compute a PD of A with orthogonal rank-1 terms 

for all nonzero Oj do 

Ay* 4— reshape into a nf 1 x • • • x tensor 

end for 
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The PD with orthogonal rank-1 terms is easily computed for the case d = 2 as 
the SVD. When d > 3, several options are available. A first option is to compute 
the CPD with the additional constraint of orthogonal factor matrices. This orthog¬ 
onality constraint limits the size of the factor matrices and consequently also the 
total number of rank-1 terms that are possible to find. As a result, the CPD with 
orthogonality constraints does not lend itself very well to applications. We demon¬ 
strate this with a worked out example in Section [7] Alternatively, one could compute 
the HOSVD or the TTrlSVD of A, as these decompositions have orthogonal rank-1 
terms. The CPD with orthogonal factor matrices and the HOSVD can be computed 


in Matlab using Tensorlab [16], freely available from 

http : //www.tensorlab.net/ 

A Matlab/Octave implementation of Algorithm 
works for any arbitrary degree d and tensor order 

4.1 

k 

that uses the TTrlSVD and 
:an be freely downloaded from 


https://github.com/kbatseli/TKPSVD 

In Section [5j we introduce a new framework in which many different structured 
tensors (symmetric, persymmetric, centrosymmetric, Toeplitz, Hankel,...) can be de¬ 
scribed and then prove that the TKPSVD algorithm guarantees to preserve these 
structures in the KP factors. But first, we present a small modification of Algorithm 
|4.1| for the case of diagonal tensors. 


4.1. Diagonal tensors. A diagonal tensor T> is an extremely simple symmetric 
tensor (see Section [5] for the definition). If we define the main diagonal of a cubical 
tensor as the entries Aj, i 2 ... lk with i\ = *2 = • • • = then the entries not on the 
main diagonal of a diagonal tensor are per definition zero. It is easy to see that the 
Kronecker product of two diagonal tensors is also diagonal. This motivates us to ad¬ 
just Algorithm 4.1 such that only the main diagonal entries 2?, 1 j 1 ...j 1 are considered. 
This reduces the number of entries to store in memory from n d to n. As a result, the 
diagonal tensor V is decomposed into a Kronecker product of diagonal factors V^' 1 . 
Suppose a degree d TKPSVD of a diagonal tensor V is required. We then consider 
the vector a that contains all main diagonal entries with entries aq^...^ and reshape 
it into a d -way tensor A. Note that since the indices are already in the right order, 
no permutation of indices is required and the PD decomposition can be directly com¬ 
puted from A. Each KP factor 1 of (1.2) is then an rq x ■ ■ ■ x rq diagonal tensor 
with main diagonal entries a - . The pseudocode for the diagonal TKPSVD algorithm 


is summarized in Algorithm 4.2 


Algorithm 4.2. TKPSVD Algorithm for a diagonal tensor 
Input: diagonal tensor D, dimensions ni, ri 2 ,..., nd 

(i) 

Output: a i,. • •, ctr , diagonal tensors D ■ 
a <r- main diagonal entries of D 
A <r- reshape a into an n\ x 712 x ... x nd tensor 

a^\ . .., ajp, <7i,..., (Jr t— compute a PD of A with orthogonal rank-1 terms 
for all nonzero Oj do 

T>j ^ 4 — a ni x • • • x ni diagonal tensor with main diagonal entries 

end for 

It is interesting to investigate whether it is possible to adjust Algorithm |4.1| to exploit 
other specific structures, like the general symmetries which we define in Section [5j At 
first sight this does not seem to be straightforward to exploit, since A will not retain 
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the original structure. We keep this problem for future research. 



5. General symmetric tensors. It turns out that the d- way tensor A in Algo¬ 
rithm [47T] allows us to generalize the notion of symmetric tensors in a very natural way. 
The motivation of introducing general symmetry lies in the fact that then only one 
proof suffices to show the preservation of symmetry, persymmetry, centrosymmetry 
and many other symmetries in the KP factors of the TKPSVD. This new framework 
also provides a different perspective of describing and investigating these different 
tensor structures. Figure 4.3 shows an overview of how the notion of general sym¬ 
metry encapsulates symmetric, persymmetric, centrosymmetric, Toeplitz and Hankel 
tensors. The key idea of the general symmetric structure is that it involves particu¬ 
lar permutations P of the entries of vec(A) that can be decomposed into Kronecker 
product of smaller permutations along each mode of A. We discuss and demonstrate 
this decomposition of the permutation P for three particular cases. In the remainder 
of this Section we always assume that A is a k way cubical tensor of dimensions n. 


5.1. Symmetry. The symmetric structure of a k- way cubical tensor A can be 
defined as a particular permutation of the entries of vec(A). This permutation is 
described by the perfect shuffle matrix. 

Definition 5.1. The perfect shuffle matrix S is the n k x n k permutation matrix 


/ 7(1 : n k ~ x : n k ,:) \ 
7(2 : n k ~ x : n k ,:) 


\I(n k ~ 1 : r ? fc_1 : n k ,:)J 


where I is the n k x n k identity matrix and Matlab colon notation is used to denote 
submatrices. 

It is easily verified that for a symmetric k- way tensor A we have that S vec(A) = 
vec(_4). We now turn this reasoning on its head and define a symmetric tensor as 
any tensor A that satisfies S'vec(A) = vec(A). The perfect shuffle matrix S reduces 


generalizes the notion of a perfect shuffle matrix to tensors. 


to the matrix defined in [§] p. 86] for the case k = 2. In this sense, Definition 5.1 
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In what follows we will apply the TKPSVD algorithm to construct the A tensor 
and see how this affects the equation Svec(./4) = vec(*4). In order to illustrate this 
process we will consider the 3-way example tensor from Section [4] and suppose that it 
is symmetric. This symmetry implies that 


(5.1) +*1*4*7] [* 2 * 5 * 8 ] [*3*6*9] — • / 4 [* 2 * 5 * 8 ] [*3*6*9] [*1*4*7] *^[*3*6*9] [* 1*4 *7] [* 2*5 * 8 ] * 


In other words, the symmetry of A is equivalent with swapping i\ with either ig or * 3 , 
»4 with either i 5 or i 6 and *7 with either i 8 or ig. The TKPSVD algorithm reshapes and 
permutes the symmetric tensor A into the tensor A, with entries [* 4 i s i e ] [± 7 £s» 9 ] * 

Although A is not symmetric, the symmetry of A still allows us to swap the indices 


as indicated in (5.1l such that 


•^•[*1*2*3] [*4*5*6] [*7*8*9] — +*4*5*6] [*7*8*9] [* 1 * 2 * 3 ] — — '4-[*7* 8 *g][*l*2*3][*4*5*6]l 

which can be rewritten as 


(5.2) 


A X! 5*1 x 2 S 2 X 3 S 3 = A 


where all 5) matrices are perfect shuffle matrices. By using (2.1), equation (5.2) can 
be rewritten as 


(5.3) (S 3 g S 2 g Si) vec(+ = vec(+, 

which is nothing else but a reformulation of the symmetry S vec(+ = vec(+ in terms 
of vec(+. If Q denotes the permutation matrix such that Qvec(A) = vec(+, then 
from 


(S 3 g £2 g Si) vec(+ = vec(+, 

<*=> (S 3 g S 2 g Si) Qvec(A) = Qvec(A), 

<=> Q t (S 3 g S 2 g Si) Q vec(A) = vec(A), 

we infer that S = Q T (S 3 g) S 2 <S> Si) Q. This can be interpreted as the perfect 
shuffle matrix S being “decomposed” into a Kronecker product of smaller perfect 
shuffle matrices. Another way of seeing this equality is that S and S 3 g S 2 g Si are 
permutation similar. 

5.2. Centrosymmetry. Another interesting and useful permutation of the en¬ 
tries of vec(A) is the exchange matrix J, which is the n k x n k column-reversed identity 
matrix. This permutation maps each index i j of vec(A) to n — ij + 1, e.g., for the 
3-way tensor A from Section 0 the entry vec(A)[i li 4 i 7 i 2 i 5 i 8 i 3 i 6 i 9 ] is mapped to 

vec(* 4 )[n—H+l n —24 + 1 n—^7+1 n — ^2+1 n—25+1 n—28 +1 n—23+I n — 26+1 n—29+I] 

and vice-versa. A k- way cubical tensor A is defined to be centrosymmetric when 

Jvec(A) = vec(A). 

Our definition of centrosymmetric tensors is equivalent with the alternative definitions 
given in j3] [T9]. The “decomposition” argument of J is completely analogous to 
the decomposition of the perfect shuffle matrix S for symmetric tensors. Following 
the TKPSVD reshapings and permutations leads to the expression J = Q T (J 3 g 
Ji g Ji) Q, where Ji, J 2 , J 3 are exchange matrices and Q is the same permutation as 
described in Subsection 15. II 
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5.3. Persymmetry. Given the definition of the perfect shuffle matrix S and 
exchange matrix J, we now define a k -way cubical tensor A to be persymmetric when 

SJvec(A) = vec(A), 

applies. Using similar arguments as in the symmetric and centrosymmetric cases, we 
can write the following decomposition S J = Q T (S 3 J 3 <g> 62 J -2 ® Si J\) Q. Using the 
mixed-product property of the Kronecker product, we can rewrite the permutation 
decomposition as S J = Q T (S 3 ® S 2 <S> <S'i)(J 3 ® Ji)Q- 

5.4. General symmetric tensor. We now define general symmetric tensors by 
generalizing the previous three examples of particular symmetries. 

Definition 5.2. A k-way cubical tensor A is a general symmetric tensor if 

P vec(A) = vec(A ), 

where the permutation matrix P can be written for any arbitrary degree d and di- 
mensions n K r J into a Kronecker product of smaller permutation matrices P \,..., Pd 
as 


(5.4) 


P = Q T (P d ®---®P 2 ®Pi)Q, 


where Q is the permutation matrix such that Q vec(A) = vec(A) and the dimension 
of each of the permutation matrices Pi is nJu^ 0 - 


General skew-symmetric tensors are defined similarly as in Definition 5.2 where 
now P vec(A) = —vec(.A) needs to hold. The permutation matrices Pi,..., Pfc do not 


necessarily need to be of the same dimension. In fact, the definition requires that (5.41 


holds for any arbitrary degree d and dimensions ni l> of the KP factors. For example, 
there are several ways in which a general symmetric 12 x 12 x 12 tensor A can be 
decomposed into a TKPSVD. There are three different decompositions when d = 3, 
depending on the order of the 3x3x3,2x2x2 and 2x2x2 tensors. Likewise when 
d = 2 there are different orderings of the 3x3x3,4x4x4 or 6x6x6, 2x2x2 tensors. 
Each of these TKPSVDs is characterized by different Pfc and Q permutation matrices, 


nevertheless, (5.4) needs to hold for all of them for A to be general symmetric. 


5.5. Shifted-index structure. Within the set of general symmetric tensors 
there are other interesting, more restrictive structures that we call shifted-index struc¬ 
tures. These are tensors whose entries do not change when at least one index is 
“shifted”. 

Definition 5.3. A k-way cubical tensor A has a shifted-index structure if 


^[ii][i2]---[ifc] ^[ii+Ai][i2+A 2 ]---[ifc+A ( .] j 

where at least one of the integer shifts Ai,..., Afc is nonzero. For any two nonzero 
shifts A i, A j, either A; = A j or A i = — Aj must be satisfied. 

Of course, none of the shifted indices i\ + A 1; ..., i k + A k are allowed to go “out 
of bounds”. The case where Ai = A 2 = • • • = Afc is called a Toeplitz tensor and is 
a special case of a persymmetric tensor. Similarly, a symmetric tensor for which all 
shifts are zero except for one arbitrary pair Aj = — Aj is called a Hankel tensor. A 
tensor for which all shifts are zero except A r has constant entries along the r fibres. It 
is straightforward to show that any shifted-index structure is also a general symmetry 
by writing down its corresponding permutation matrix and showing that Definition 
|5.2| applies. 
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(b) Hankel permutation matrix P 


Fig. 5.1: Permutation matrices for a 27 x 27 x 27 Hankel tensor. 


Example 5.1. Consider a 27 x 27 x 27 Hankel tensor H. The Hankel structure 
means that the tensor is also symmetric, which implies that S vec(TL) = vec(H) with 
S the 19683 x 19683 perfect shuffle matrix. Now consider a degree-3 TKPSVD where 
each KP factor is a 3x3x3 tensor. We can then retrieve S from the 27 x 27 perfect 
shuffle matrix Si = S 2 = S 3 as S = Q T (S 3 ® S 2 <8> Si) Q, where Q is the permutation 
matrix in vecflH) = Q vec(H). The 27 x 27 permutation matrices P\ = P 2 = P 3 that 
define a 3 x 3 x 3 Hankel tensor A are completely specified by the vector of indices 


i = [1,4,5,10, 7,8,11,12,15,2,13,14,19,16,17,20,21,24,3,22,23,6,25,26,9,18,27], 


since vec(A)(i ) = vec(A), where vec(A)(i) is Matlab notation to denote P 3 vec(A). 
If we now set P = Q T (P 3 ® P% ® Pi) Q then indeed P vec(H) = vec(TL) is satisfied. 
Figure \5.1\ shows the nonzero pattern of both S and P. Observe that although the 
Hankel permutation P is a special case of a symmetry, the nonzero pattern is very 
different from that of S. 


6. Preservation of structures. It is quite a remarkable fact that all general 
symmetries, included the shifted-index structures, are prese rved in the cubical KP 

The orthogonality 


factors ^4^ 


4.1 


when they are computed according to Algorithm 
of the rank-1 terms in the PD plays a crucial role in this. Another critical element 
are the scalar coefficients crj’s of the TKPSVD, which are required to be distinct. We 
now show how this comes about. 


6.1. General symmetry. In order to prove general symmetry preservation in 
the KP factors, we need the following useful lemma. 

Lemma 6.1. Suppose a = vec(A) £ R n xl with a T a = 1 and P is a permutation 
matrix that corresponds with a general symmetry. Then the cubical tensor A obtained 
from reshaping a is general symmetric if and only if a T P a = 1 or general skew- 
symmetric if and only if a T P a = —1. 

Proof. We first prove Pa = a => a T P a = 1. Since a has unit norm we can 
write a T a = 1 and substitution of a by Pa then results in a T Pa = 1. The proof for 
a T P a = 1 =£* Pa = a goes as follows. Let b = P a, then we have that || 6||2 = 1 and 
a T b = cos a. Since cos a = 1 it follows that a = 0 and a is a multiple of b , but since 
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||a|| 2 = H&H 2 = 1 it follows that a = b = P a. The proof for the skew-symmetry of A 
follows the same logic. □ 

The general symmetry of A implies that 


( 6 . 1 ) 


A = A x 1 Pi x 2 Pi x 3 • • ■ x d P d , 


where all Pf s are permutation matrices. We now substitute A in both sides of (6.11 
by its PD and obtain 


( 6 . 2 ) 


R 

E (l) (2) (d) 

cr J ay O a) O ■ ■ ■ o a,j = 

l=i 


R 

i— 1 


Piaf ] o P 2 aj 2) o • 


■ • ° Pda 


(d) 


The orthogonality of each rank-1 term and | 2 = 1 implies that the mode products 

of both sides of ( |6.2| ) with {a ^) T ,..., (a^) T along the modes 1,2 ,d respectively 
for any ft e {1,...,/?} C N results in 


(6.3) 




= X>j (4 ,, ) T -r .<’ K’T'PiaY' ••• 

J=1 


J 2 ) 




(d) 


We have that each of the (a^ ) T Pia^ 

,(0i 


scalars lies in the real interval [—1,1], since 


i \ay 11 2 = 1 for all j and P is a permutation matrix. We now assume that the following 
condition holds. 

Condition 1. All aj's in the PD of A are distinct and all terms on the right-hand 
side of (6.3) except for the one corresponding with <Jk vanish. 

The equality in (6.3) holds under Condition [l] when (o^) T Pi a^' 1 = 0 for at least 
one particular i when k ^ j and when 


(6.4) 


n ^kfp^k = i- 


The only possible way for (6.4) to be true under the constraint that each of the 
( a k) T Pi a ^k sca l ars H es in the interval [—1,1] is when 


(6.5) 


KTPi 




= ± 1 . 


From Lemma 


6.1 


we know that if the right-hand side of (6.5) is 1, then A[^ is general 
symmetric, otherwise A^ is general skew-symmetric. In addtion, (6.4) implies that 
there are either zero or an even number of general skew-symmetric KP factors in 
each term. This proves the main theorem on general symmetry-preservation in the 
TKPSVD. 


Theorem 6.2. Let A be a general symmetric tensor with a dth-degree TKPSVD 
into cubical KP factors If Condition [l] holds, then each of the factors in the 

TKPSVD is either a general symmetric or a general skew-symmetric tensor. There 
are always either zero or an even number of skew-symmetric factors in each term of 

Let us see what happens when Condition [1] is not satisfied. Suppose that cr^ = 
<Jk+i and that the terms corresponding with the other <jj’s vanish. We can then 


( 1 . 2 ). 
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combine the <Jk and ctk+i terms such that now 


( 6 . 6 ) 




= i 


needs to hold. Contrary to the case of distinct ctj’s, there are now multiple ways 


that (|6.6[) can be satisfied without {a'i > ) T Pi a^ = ±1 being true, which implies that 


the terms corresponding with ak,<Jk+i will not necessarily have general symmetric 
factors. A very particular case where Condition [I] is not satisfied for all terms is for 
symmetric tensors of order k > 2. We discuss this case, along with other general 
symmetries, in Section [7] 

7. Numerical Experiments. In this section we discuss numerical experiments 
that illustrate different aspects of the TKPSVD algorithm. We will demonstrate the 
structure preservation of generaly symmetries and shifted-index structures and discuss 
a curious observation for symmetric tensors. We also compare the use of the CPD 
with orthogonal matrix factors, the HOSVD and TTrlSVD in terms of runtime and 
storage. Finally, we illustrate how the Kronecker product structure can be interpreted 
as a multiresolution decomposition of images. All computations were done in Matlab 
on a 64-bit 4-core 3.3 GHz desktop computer with 16 GB RAM. 

7.1. General symmetric structure. As a first example, we demonstrate the 
use of the CPD with orthogonal factor matrices, the HOSVD and the TTrlSVD to 
compute the TKPSVD, together with the preservation of centrosymmetry in the KP 
factors. 

Example 7.1. We construct a 24 x 24 x 24 centrosymmetric tensor A with its 
distinct entries drawn from a standard normal distribution and compute a 3rd degree 
decomposition with factor sizes 4x4x4,3x3x3, 2x2x2 respe ctive ly. The reshaping 


and permutation steps result in a 8 x 27 x 64 tensor A. Table 7.1 compares the use 


of the CPD with orthogonal factor matrices, the HOSVD and TTrlSVD for the PD 
of A. We list the total number of rank-1 terms, the required memory for storage of 
the PD of A, the total runtime to compute the TKPSVD (computed as the median 


over 100 runs) and the relative error ||A — X)p=i total 

number of rank-1 terms for the orthogonal CPD is limited to only 8, since the first 
factor matrix will be an orthogonal 8x8 matrix. This results in a large relative error. 
The runtime for computing the orthogonal CPD is also very long compared to the 
HOSVD and TTrlSVD. The main difference between the HOSVD and TTrlSVD lies 
in the total number of rank-1 terms. But since the HOSVD reuses the mode vectors, 
this results in slightly less required memory. All rank-1 terms computed with both the 
HOSVD and TTrlSVD retain the centroysymmetric structure. The TTrlSVD method 
results in 56 terms that have 2 skew-centrosymmetric factors, The HOSVD has 1792 
such terms. Note that the HOSVD has a 8 x 27 x 64 core tensor containing 6912 
nonzero entries. 

Due to the limitations of the CPD with orthogonal factors, as demonstrated in 
Example |7.1[ we will refrain from using it for the following examples in this Section. 
We now demonstrate the occurrence of cr’s with multiplicities for symmetric tensors. 
Consequently, KP terms that belong to the same o will not inherit the general sym¬ 
metry in their factors. This is true for when both the TTrlSVD and HOSVD are 
used. The TTrlSVD case however has a lot more regularity than the HOSVD. When 
the TKPSVD of a k -way symmetric tensor is computed with the TTrlSVD, then all 
multiple cr’s have a multiplicity of k — 1. 


id) i 
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Table 7.1: Comparison of CPD, HOSVD and TTrlSVD for the TKPSVD. 


Method # Terms Storage (kB) Runtime (seconds) Relative Error 


orthogonal CPD 8 3.14 

HOSVD 6912 146.19 

TTrlSVD 216 154.06 


768.07 

0.18 

0.21 


0.947 

2.21e-15 

2.39e-15 


Example 7.2. Consider a symmetric 8x8x8 tensor with distinct entries drawn 
from a standard normal distribution. We compute its degree-3 TKPSVD using the 
TTrlSVD and obtain 56 KP terms. For this 3-way tensor, each multiple a has a 
multiplicity of k — 1 = 2. There are 8 such pairs, which implies that 16 terms are not 
(skew)-symmetric. Next, we compute a degree-3 decomposition for an 8x8 symmetric 
matrix, which results in 14 KP terms. All the a’s are distinct, which implies that all 
terms in the decomposition are (skew)-symmetric. Finally, we compute the degree-3 
TKPSVD of a 4-way symmetric tensor. This decomposition consists of 230 terms. 
There are 20 3-tuples of multiple a’s, which means that 60 KP terms are not (skew)- 
symmetric. 

7.2. Shifted-index structure. Next, we investigate the dependence of the total 
number of KP terms and total runtime on the ordering of the KP factors for both the 
TTrlSVD and HOSVD. 

Example 7.3. Consider a 64 x 64 x 64 x 64 Hankel tensor with its distinct 
entries drawn from a standard normal distribution. We compute its TKPSVD into 3 
Kronecker factors with dimensions 2x2x2x2,4x4x4x4 and 8 x 8 x 8 x 8 over 
all possible orderings of the factors and investigate the total number of obtained KP 
terms for both the TTrlSVD and HOSVD, together with total runtimes. The results 
are shown in Table | 7. 2\ The first thing to notice is that the total number of terms in 
the TKPSVD and total runtime are quite independent from the factor ordering when 
the HOSVD is used. On average, about 2000 terms are needed and the computation 
takes a little over 4 minutes. The TTrlSVD needs about 20 times less terms and is 
for one particular ordering more than 80 times faster. Note that although the HOSVD 
requires more terms, just like in Example | 7. 1\ it will require less memory for storage 
due to the fact that it reuses the mode vectors for each KP term. All of the terms in 
every of the decompositions retained the Hankel structure. 

Table 7.2: Number of KP terms and total runtime for TTrlSVD and HOSVD. 


Ordering 

ff Terms 

Runtime (seconds) 


TTrlSVD 

HOSVD 

TTrlSVD 

HOSVD 

2,4,8 

65 

1968 

9.76 

254.30 

2,8,4 

65 

1973 

2.47 

251.70 

4,2,8 

65 

1993 

5.80 

254.96 

4,8,2 

65 

2146 

5.75 

256.41 

8,2,4 

145 

2067 

247.66 

249.97 

8,4,2 

145 

2105 

239.95 

254.16 
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Fig. 7.1: Original 4000 x 6000 x 3 image. 


7.3. Multiresolution decomposition of images. An interesting illustration 
of the TKPSVD is in the multiresolution decomposition of a m x ri 2 x 3 colour image 
A. This example gives an interpretation to two different aspects of the TKPSVD: 
truncation of the Kronecker product and truncation of the number of terms. Indeed, 


every pixel of the n[ d> x x 3 image A^ d> is “blown up” by each Kronecker product 
in A ® A^ d ~^ ® • • • <S> A 1 - 1 ' 1 until the resolution ni x ri 2 x 3 is obtained. Truncating 
the Kronecker products to only a few factors hence effectively reduces the resolution. 
Furthermore, compression can be achieved at different resolutions by retaining only 
a few terms. The compression rate achieved for retaining k factors and r number of 
terms in (1.21 is defined as 




k— 1 ( d—i) ( d—i) (d— 1 


n 


rin 


Tl O 


k— 1 (d—i) (d—i) (d— 1) ' 


E k — 

i=i 


n. 


This is illustrated with the 6000 x 4000 x 3 colour image in Figure [7~i[f~| A degree-5 
TKPSVD is computed with dimensions 


(250 x 375 x 3) ® (2 x 2 x 1) ® (2 x 2 x 1) ® (2 x 2 x 1) ® (2 x 2 x 1). 


The total runtime to compute the TKPSVD using the TTrlSVD and HOSVD was 18 
and 30 seconds respectively. In contrast, computing a standard SVD of one of the three 
slices A (:,:, i) takes about 1 minute and consists of 4000 rank-1 terms. The TTrlSVD 
needs 240 terms while the HOSVD needs 65536. For this particular example, the 
compression rate when retaining k factors and r terms can be approximated by 

(250 x 375 x 3) (2 x 2 x l) fe " 1 _ (2 x 2 x l) fe_1 _ 4 fc " 1 
r(250 x 375 x 3 + (k - 1)(2 x 2 x 1)) ~ r “ r ' 

This implies that the maximal compression rate, when r = 1, is dependent on the 
resolution, viz. the number of KP factors k in each term. At the largest resolution 


1 absfreepic.com/free-photos/download/water-nature-fall-6000x4000_90673.html 
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{k = 5), the maximal compression rate is approximately 256 while at the smallest 
resolution no compression is possible through truncation of KP terms. A common 
measure to quantify the quality of reconstruction of lossy compressed images is the 
peak signal-to-noise ratio (PSNR). The PSNR is defined as 

PSNR = 201og 10 (MAXi) - 101og 10 (MSE) 


where MAXi is the maximal possible pixel va lue, 255 in our case, and MSE is the mean 
squared error ||A — A\\ 2 F /{n\ ■ • 3). Figure 


7.2 


shows the PSNR as a function of the 
number of retained KP terms, computed from the TTrlSVD for the highest possible 
resolution. In this case, the PSNR can be completely determined from the a’s in the 
TKPSVD using ( |1.4| ). Acceptable values of the PSNR are between 30 and 50 dB and 
are obtained from retaining the first 20 terms, which corresponds with a compression 
rate of about 12.8. Figure [773] displays 1-term approximants for 3 different resolutions. 
For Figures [A3t a),(b),(c) the PSNR is 58dB, 58dB and 57dB respectively. 



Fig. 7.2: Increase of the PSNR as the number of TTrlSVD-KP terms grows for the 
4000 x 6000 x 3 resolution. 


8. Conclusions. In this paper, we introduced the tensor Kronecker product 
singular value decomposition that decomposes a real k -way tensor A into a linear 
combination of tensor Kronecker product terms with an arbitrary number of d fac¬ 
tors A = a j -A 1 / 1 ® ® ~^P ■ This decomposition enables easy computation 

of a Kronecker product approximation and a very straightforward determination of 
the relative approximation error without explicit construction of the approximant. 
We proved that for many different structured tensors, the Kronecker product factors 
Ap ,..., Ap are guaranteed to inherit this structure. In addition, we introduced the 
new framework of general symmetric tensors, which includes many different structures 
such as symmetric, persymmetric, centrosymmetric, Toeplitz and Hankel tensors. 

Acknowledgements. The authors would like to thank Martijn Bousse and Nico 
Vervliet for their invaluable help on computing a CPD with orthogonal factor matrices 
in Tensorlab. 
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(a) 250 x 375 


(b) 500 x 750 



(c) 1000 x 1500 


Fig. 7.3: First term of the KP decomposition for 3 different resolutions. 


REFERENCES 

[1] K. Batselier, H. Liu, and N. Wong, A constructive algorithm for decomposing a tensor 

into a finite sum of orthonormal rank-1 terms, SIAM J. Matrix Anal. Appl., 36 (2015), 
pp. 1315-1337. 

[2] J.D. Carroll and J.-J. Chang, Analysis of individual differences in multidimensional scaling 

via an n-way generalization of “Echart-Young” decomposition, Psychometrika, 35 (1970), 
pp. 283-319. 

[3] H. Chen, Z. Chen, and L. Qi, Centro symmetric, skew centro symmetric and centro symmetric 

cauchy tensors, ArXiv e-print 1406.7409, (2014). 

[4] L. De Lathauwer, B. De Moor, and J. Vandewalle, A multilinear singular value decom¬ 

position, SIAM J. Matrix Anal. Appl., 21 (2000), pp. 1253-1278. 

[5] W. HACKBUSCH, B. N. Khoromskij, and E. E. Tyrtyshnikov, Hierarchical Kronecker tensor- 

product approximations, J. Numer. Math., 13 (2005), p. 119156. 

[6] R. A. Harshman, Foundations of the PARAFAC procedure: Models and conditions for an 

“explanatory” multi-modal factor analysis, UCLA Working Papers in Phonetics, 16 (1970), 
p. 84. 

[7] J. Kamm and J. G. Nagy, Optimal Kronecker Product Approximation of Block Toeplitz Ma¬ 

trices, SIAM J. Matrix Anal. Appl., 22 (2000), pp. 155-172. 

[8] T.G. Kolda and B.W. Bader, Tensor decompositions and applications, SIAM Rev., 51 (2009), 

pp. 455-500. 



TKPSVD 


21 


[9] C. F.Van Loan, The ubiquitous Kronecker product, J Comput Appl Math, 123 (2000), pp. 85 
- 100 . 

[10] C. F. Van Loan and N. Pitsianis, Approximation with Kronecker Products, in Linear Algebra 

for Large Scale and Real Time Applications, Kluwer Publications, 1993, pp. 293-314. 

[11] J.G. Nagy and M.E. Kilmer, Kronecker product approximation for preconditioning in three- 

dimensional imaging applications, IEEE T. Image Process., 15 (2006), pp. 604-613. 

[12] J. G. Nagy, M. K. Ng, and L. Perrone, Kronecker Product Approximations for Image 

Restoration with Reflexive Boundary Conditions, SIAM J. Matrix Anal. Appl., 25 (2003), 
pp. 829-841. 

[13] P. A. Regalia and M. K. San jit, Kronecker Products, Unitary Matrices and Signal Processing 

Applications, SIAM Review, 31 (1989), pp. 586-613. 

[14] M. Rezghi, S. M. Hosseini, and L. Elden, Best Kronecker Product Approximation of The 

Blurring Operator in Three Dimensional Image Restoration Problems, SIAM J. Matrix 
Anal. Appl., 35 (2014), pp. 1086-1104. 

[15] J. Salmi, A. Richter, and V. Koivunen, Sequential unfolding svd for low rank orthogonal 

tensor approximation, in Signals, Systems and Computers, 2008 42nd Asilomar Conference 
on, Oct 2008, pp. 1713-1717. 

[16] L. Sorber, M. Van Barel, and L. de Lathauwer, Tensorlab Version 2.0. Available online, 

January 2014. 

[17] L. Tucker, Some mathematical notes on three-mode factor analysis, Psychometrika, 31 (1966), 

pp. 279-311. 

[18] L. R. Tucker, The extension of factor analysis to three-dimensional matrices, in Contributions 

to mathematical psychology., H. Gulliksen and N. Frederiksen, eds., Holt, Rinehart and 
Winston, New York, 1964, pp. 110-127. 

[19] X. Zhao and Q. Yang, The spectral radius of nonnegative centrosymmetric tensor, J. High 

School Numer. Math., 36 (2014), pp. 58-66. 



